Close-packed floating clusters: granular hydrodynamics beyond the freezing point? 
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Monodisperse granular flows often develop regions with hexagonal close packing of particles. 
We investigate this effect in a system of inelastic hard spheres driven from below by a "thermal" 
plate. Molecular dynamics simulations show, in a wide range of parameters, a close-packed cluster 
supported by a low-density region. Surprisingly, the steady-state density profile, including the close- 
packed cluster part, is well described by a variant of Navier-Stokes granular hydrodynamics (NSGH) . 
We suggest a simple explanation for the success of NSGH beyond the freezing point. 
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Continuum modeling of flow of macroscopic grains re- 
mains a challenge The best known version 
of continuum theory here is the Navier-Stokes granular 
hydrodynamics (NSGH) for a system of inelastic hard 
spheres @, 0- The applicability of NSGH is limited 
to rapid granular flows By definition, these flows 
are dominated by binary particle collisions, while multi- 
particle interactions are negligible. Despite this drastic 
simplification, the validity of the NSGH demands sev- 
eral additional assumptions, some of which can be rather 
stringent. Under the molecular chaos assumption, the 
NSGH is derivable systematically from more fundamen- 
tal kinetic equations for inelastic hard spheres 0, 0- E} • 
Going over from kinetic equations to hydrodynamics, one 
should assume scale separation: the mean free path of the 
particles should be much smaller than the characteristic 
length scale, and the mean time between two consecu- 
tive collisions much shorter than any characteristic time 
scale, described hydrodynamically. The inelasticity of 
particle collisions brings immediate complications. Al- 
ready at moderate inelasticity q = (1 — r)/2 (where r is 
the coefficient of normal restitution of the particle colli- 
sions) , the scale separation may break down, even in the 
low-density limit [111 Il2| . The normal stress difference 
12] and deviations of the particle velocity distribution 
llL Hii ] from the Maxwell distribution also become im- 
portant for moderately inelastic collisions. Therefore, the 
NSGH is expected to be accurate only for small inelas- 
ticity, q <C 1. 

Additional complications appear at large densities. 
Here the molecular chaos assumption breaks down, al- 
ready for elastic hard spheres, when the packing frac- 
tion approaches the freezing point value 4>f — 0.49 (in 
three dimensions) or 0.69 (in two dimensions). As the ki- 
netic equations become invalid, the constitutive relations 
(CRs), necessary for the closure of hydrodynamics, are 
not derivable from first principles anymore. This is the 
regime considered in this work. We consider an ensemble 
of monodisperse, nearly elastic hard spheres in such con- 
ditions that the standard NSGH @,uj breaks down be- 
cause of large densities, not large inelasticity. Our main 



objective is to check whether a variant of NSGH can still 
be used in an extreme case when the packing fraction is 
close to the maximum possible value, corresponding to 
hexagonal packing of spheres. 

We will focus on granular materials fluidized by 
a rapidly vibrating bottom plate in a gravity field. 
Vibrofluidized granular materials exhibit fascinating 
pattern-formation phenomena that have attracted much 
recent interest 14]. In the high-frequency and small- 
amplitude limit of vibrofluidization, there is no direct 
coupling between the vibration and the collective granu- 
lar motion. In a simplified description of this limit one 
specifies a constant granular temperature at an immobile 
bottom plate. In a wide range of parameters, molecu- 
lar dynamics (MD) simulations of this system show an 
(almost) close-packed cluster of particles, floating on a 
low-density fluid, see below. The close-packed floating 
cluster is an extreme form of the density inversion, a 
phenomenon well known in vibrofluidized granular ma- 
terials. Lan and Rosato [w| were apparently the first to 
observe density inversion in three-dimensional MD sim- 
ulations. Kudrolli et al. 0] observed a floating cluster 
in a reduced-gravity experiment: a slightly tilted two di- 
mensional system of steel spheres rolling on a smooth 
surface and driven by a vibrating side wall. Recently, a 
pronounced density inversion has been observed, in two- 
and three dimensional vibrofluidized granular beds, by 
Wildman et al. [l7j . 

An accurate hydrodynamic description of almost close- 
packed floating clusters seems a very difficult task, as the 
packing fraction here is far beyond the freezing point. 
Still, we will attempt to use a variant of NSGH for this 
purpose. This attempt will prove successful, and we will 
suggest an explanation. Here is the model problem we 
are working with. Let N 3> 1 nearly elastic hard spheres 
of diameter d and mass m move in a two-dimensional box 
with periodic boundary conditions in x-direction (period 
L x ) and infinite height. The driving base is located at 
y = 0. Gravity acceleration g acts in negative y direc- 
tion. Upon collision with the base, the particle velocity is 
drawn from a Maxwell distribution with temperature Tq 
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(which is measured in the units of energy). The kinetic 
energy of the particles is being lost by inelastic hard- 
core collisions parameterized by a constant inelasticity 
parameter q -C 1. Figure shows a typical snapshot of 
an almost close-packed floating cluster, observed in an 
event-driven MD simulation of this system. Hexagonal 
packing is apparent in Fig. HI (lif. 
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FIG. 1: A snapshot showing the close-packed floating cluster. 
The parameters are N = 10 4 , L x = 100, T = 2\/3 ■ 10 2 , 
r — 0.98815 and g = d = m = 1. The figures on the right are 
magnifications of the indicated areas. 



Going over to a hydrodynamic description of zero- 
mean-fiow steady states, we introduce coarse-grained 
fields: the particle number density n(y), the granular 
temperature T(y) and the pressure p(y). The maximum 
possible value of n is the hexagonal close-packing value 
n c = 2/(v / 3d 2 )- A laterally uniform steady state is de- 
scribed by the momentum and energy balance equations: 
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where k is the thermal conductivity and / is the energy 
loss rate by collisions. To proceed, we need CRs: an 
equation of state (EOS) p = p(n, T) and relations for 
k and / in terms of n and T . First-principles CRs are 
available only in the low-density limit (well below the 
freezing point). Grossman et al. derived a set of 
approximate global CRs for a version of NSGH that as- 
sumes nearly elastic collisions, but is not limited to low 
densities. Grossman et al. employed free volume argu- 
ments in the close vicinity of the hexagonal packing and 
suggested simple interpolations between the hexagonal- 
packing limit and low-density limit. These interpolations 
include two fitting constants a and 7 (see below). The 
optimum values of these constants were found by a com- 
parison with MD simulations of a system of inelastic hard 
spheres driven by a thermal wall at zero gravity [Tlj . 

Notice that, prescribing global CRs of any type, one 
grossly simplifies the delicate issue of phase coexistence 



that is expected to occur here in close analogy to the 
system of elastic hard spheres 0, • Still, we will use 
the simple CRs to attempt a NSGH description of 
the close-packed floating clusters. In our notation, the 
CRs E3 read 
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and I = 4(^/70 gnm^T 3 / 2 . 
path of the grains, 

I = 
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Here I is the mean free 
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and a = 1 - (3/8) 1 ' 2 
a = 1.15 and 7 = 2.26 



According to Grossman et al. 
We adopted this value of 7, but 
found better agreement between the hydrodynamics and 
MD (see below) for a = 0.6. The value of fjt = 0(1) 
drops out from the steady-state problem. 

Recently a more accurate global EOS p = p(n, T) has 
been suggested |2cij. Still, in the absence of comparably 
accurate relations for k and /, employing a more accurate 
EOS in Eqs. would be an excess of accuracy. 

Equations should be complemented by three 
boundary conditions. One of them is T(0) = Tq = const. 
Integrating the first of Eqs. JU over the height from 
to 00 and using the conservation of the total number of 
particles: J„ n(y) dy = N/L x = const, we obtain the 
second boundary condition: p(0) = mgN/L x . The third 
one is a zero heat flux (that is, a constant granular tem- 
perature) at y — > 00 [2lL I22I] . In practice, one should use 
the shooting method, varying the heat flux —ndT/dy at 
y = until the third condition is satisfied with desired 
accuracy. 

Let us measure y in units of the gravity length scale 
A = To/(mg) (note that X/d should be large enough to 
ffuidize the granulate at the bottom). We rewrite Eqs. 
ijljl. in scaled form, as three first-order equations: 
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Here Z(y) = n c /n(y) is the inverse scaled density, P(y) = 
p(y)/(n c To) is the scaled pressure, and — $(y) is the 
scaled heat flux. The functions F\ , F2, and Q are 



F 1 (Z) = 



aZ(Z - 1) + ^/32/3(Z - a) 



(Z - 1) 3 / 2 Z 3 / 2 
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Finally, A = (64/ 7 ) q (X/d) 2 and / = {VSd 2 N)/(2 XL X ) 
are two scaled governing parameters. Parameter A con- 
trols the relative role of the inelastic heat losses and heat 
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FIG. 2: The critical value A c of the hydrodynamic inelasticity 
parameter A for a density inversion versus the effective area 
fraction /. At small / A c scales like f~ 2 . 



conduction, while / is the effective area fraction of the 
grains (it can be smaller or greater than unity). The 
boundary conditions at the base become 



P(0) = / and Z(0) 



l + / + (l + 6/ + / 2 ) 1/2 
2/ 



(6) 



Using the hydrodynamic formulation, we first deter- 
mine the condition for a density inversion. At too small 
inelasticity q (the rest of the parameters fixed) there is 
no density inversion, like in the elastic case q — where 
T(y) = const and n(y) goes down monotonically. At 
large enough q the temperature T(y) drops rapidly with 
y. To maintain the hydrostatic balance, n(y) should in- 
crease with the height, on an interval of heights between 
y = and the location of the density maximum y = y c . 
In our hydrodynamic formulation the density inversion 
occurs, at fixed /, when A > A c , where A c = A c (/) is 
a critical value. The density inversion is born at y = 0: 
A = A c (/) corresponds to dn/dy vanishing at y = 0. Us- 
ing this condition together with Eqs. J5J and ©, we 
obtain 
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For a given /, Eq. prescribes the heat flux at the 
base y = that corresponds to the birth of the density 
inversion. Using shooting, we determine, for every /, the 
critical value A c (/), demanding that the temperature ap- 
proaches a constant value at large heights. This proce- 
dure yields the critical curve A = A c (/) shown in Fig. 
|5] The density inversion occurs above the critical curve 
A = A c (/), and it is more and more pronounced, at fixed 
/, as A grows. Figure (a) shows the density profiles 
at / = 0.25 and three different values of A > A c . One 
can see that, at large enough A, a hexagonally-packed 
cluster appears. The scaled parameters A = 20,015 and 
/ = 0.25 correspond to the snapshot shown in Fig. ^ 
Noticeable is a steep (exponential) density fall at the up- 
per boundary of the cluster; the exponent corresponds to 
the very low temperature there. Figure 01(b) shows the 
scaled temperature for these three cases. 
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FIG. 3: Scaled density (a) and temperature (b) versus the 
scaled height y for / = 0.25 and A = 500 (solid lines), 2000 
(dashed lines) and 20,015 (dotted lines). The dotted lines 
correspond to the snapshot shown in Fig. 
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FIG. 4: Scaled density versus the height y as predicted by hy- 
drodynamics (solid curves) and observed in MD simulations, 
for three different values of the temperature at the base. The 
parameters are described in the text. The height is measured 
here in units of d. 



Figure 01 compares the density profiles, predicted by 
this hydrodynamics (solid curves), with the profiles found 
in MD simulations with N = 10 4 particles of diameter 
d = 1, mass m = 1 and r — 0.98815. The (periodic) 
box width is L x = 100, the gravity acceleration g = 1. 
The MD simulations were done for three different values 
of the temperature at the base: To = 100\/3, 200V3 and 
300v^- The hydrodynamic parameters in these three 
cases are / = 0.5 and A = 5004; / = 0.25 and A = 20015, 
and / = 0.167 and A = 45036, respectively. One can 
see that the agreement between hydrodynamics and MD 
simulations is surprisingly good. 

An additional argument in favor of hydrodynamics fol- 
lows from the dimensional analysis of the problem. The 
full set of parameters includes d,m,r, 3, To, N and L x . 
One can always choose d = m = g — 1, so there are actu- 
ally four independent parameters. This number reduces 
to three for an x-independent steady state, as N and L x 
enter the problem only through N/L x . It is crucial that 
hydrodynamics further reduces the number of parame- 
ters: now only two scaled parameters A and / appear. 
This prediction is very robust, as it is independent of the 
particular form of the functions Fi , F2 and Q (and of the 
values of a and 7). We verified this prediction in MD sim- 
ulations by varying N, To and r, but keeping A = 20, 015 
and / = 0.25 constant. After rescaling the coordinate 
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FIG. 5: Checking the hydrodynamic scaling. The solid line 
corresponds to the MD simulation shown in Fig. 1. The 
dashed line corresponds to another simulation: with N = 
2 ■ 10 4 , L x = 100, To = 4^3 ■ 10 2 , r = 0.997051 and g = d = 
m — 1. For the dashed line the height y is shrunk by a factor 
of 2. 

y by A = To /(nig), the resulting density profiles almost 
coincide with each other, see Fig. [S] A small shift be- 
tween the two profiles, observed in Fig. |5j is apparently 
caused by small vertical oscillations of the granulate. In 
this example the oscillation amplitude is about 3 particle 
diameters 0. 

As already mentioned, the global CRs completely ig- 
nore the issue of coexistence, beyond the freezing point, 
of different phases of the granulate: the liquid-like phase, 
the random close-packed phase etc. So why are they so 
successful? We believe the reason is the following. The 
vibrofluidized steady state, considered in this work, has 
a zero mean flow. Therefore, the viscosity terms in the 
hydrodynamic equations vanish. This fact is not merely 
a technical simplification. The shear viscosity of granu- 
lar flow is finite in the liquid-like phase, and infinite in 
the (multiple) domains of the random close-packed phase. 
The effective total viscosity of the system is expected to 
diverge when the coarse-grained density slightly exceeds 
the freezing density. This invalidates any NSGH for suf- 
ficiently dense flows, and necessitates the introduction of 
an order parameter and a different type of the stress- 
strain relation into the theory, cf. Ref. Q. Luckily, 
these complications do not appear for a zero-mean-flow 
state. Indeed, the EOS, heat conductivity and inelastic 
heat loss rate do not exhibit any singularity around the 
freezing point, and all the way to the hexagonal close 
packing. Therefore, the NSGH remains reasonably accu- 
rate far beyond the freezing point. A future work should 
address the important question about the range of ap- 
plicability of the NSGH (actually, of any binary collision 
model) for solid yet vibrated phase, versus the granular 
statics approach. 
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